function [GlobalParams,gbestfitness,gbesthistory] = COA(popsize, D, xmax,xmin,vmax,vmin,nfevalMAX,func,fid)
%[gbestX,gbestfitness,gbesthistory]= IDA(PopSize,D,xmax,xmin,vmax,vmin ,MaxIter,Func,FuncId);


FEs=0;
MaxFEs=10000*D;
FOBJ = func;
i=1;
%% Optimization problem variables
lu = [xmin*ones(1,D);xmax*ones(1,D)];
VarMin      = lu(1,:);
VarMax      = lu(2,:);

%% Algorithm parameters
gbestfitness=inf;
n_coy = 5;
n_packs =8;

% Probability of leaving a pack
p_leave     = 0.005*n_coy^2;
Ps          = 1/D;
gbesthistory = nan(1, MaxFEs);
%% Packs initialization (Eq. 2)
pop_total   = n_packs*n_coy;
%costs       = zeros(pop_total,1);
coyotes     = repmat(VarMin,pop_total,1) + rand(pop_total,D).*(repmat(VarMax,pop_total,1) - repmat(VarMin,pop_total,1));
ages        = zeros(pop_total,1);
packs       = reshape(randperm(pop_total),n_packs,[]);  %分种群
coypack     = repmat(n_coy,n_packs,1);

%% Evaluate coyotes adaptation (Eq. 3)
for c=1:pop_total
    costs(c,1) = FOBJ(coyotes(c,:)',fid);
    FEs=FEs+1;
    if gbestfitness>costs(c,1)
        gbestfitness=costs(c,1);
        gbestX=coyotes(c,:);
    end
    gbesthistory(FEs)=gbestfitness;
    fprintf("COA 第%d次评价，最佳适应度 = %e\n",FEs,gbestfitness);
end
nfeval = pop_total;

%% Output variables
[gbestfitness,ibest]   = min(costs);
GlobalParams        = coyotes(ibest,:);

%% Main loop
year=0;
%while nfeval<nfevalMAX % Stopping criteria
while  1
    %% Update the years counter
    year = year + 1;
    
    %% Execute the operations inside each pack
    for p=1:n_packs
        % Get the coyotes that belong to each pack
        coyotes_aux = coyotes(packs(p,:),:);
        costs_aux   = costs(packs(p,:),:);
        ages_aux    = ages(packs(p,:),1);
        n_coy_aux   = coypack(p,1);
        
        % Detect alphas according to the costs (Eq. 5)
        [costs_aux,inds] = sort(costs_aux,'ascend');%对适应度进行排序
        coyotes_aux      = coyotes_aux(inds,:);
        ages_aux         = ages_aux(inds,:);
        c_alpha          = coyotes_aux(1,:);  %第p个种群中最好的个体
        
        % Compute the social tendency of the pack (Eq. 6)
        tendency         = median(coyotes_aux,1); %使用中位数代表种群趋势
        
        % Update coyotes' social condition
        new_coyotes      = zeros(n_coy_aux,D);
        for c=1:n_coy_aux
            rc1 = c;
            while rc1==c
                rc1 = randi(n_coy_aux);  %随机选取不等于c的随机个体
            end
            rc2 = c;
            while rc2==c || rc2 == rc1
                rc2 = randi(n_coy_aux);  %随机选择不等于rc1和c的个体
            end
            
            % Try to update the social condition according to the alpha and
            % the pack tendency (Eq. 12)
            new_c = coyotes_aux(c,:) + rand*(c_alpha - coyotes_aux(rc1,:))+ ...
                rand*(tendency  - coyotes_aux(rc2,:));
            
            % Keep the coyotes in the search space (optimization problem
            % constraint)
            new_coyotes(c,:) = min(max(new_c,VarMin),VarMax); %越界检测
            
            % Evaluate the new social condition (Eq. 13)
            new_cost = FOBJ(new_coyotes(c,:)',fid);
            
            
            
            
            
            
            % Adaptation (Eq. 14)
            if new_cost < costs_aux(c,1)
                costs_aux(c,1)      = new_cost;
                coyotes_aux(c,:)    = new_coyotes(c,:);
            end
            FEs=FEs+1;
            if gbestfitness> costs_aux(c,1)
                gbestfitness= costs_aux(c,1) ;
                gbestX=   coyotes_aux(c,:);
            end
            
            gbesthistory(FEs)=gbestfitness;
            fprintf("COA 第%d次评价，最佳适应度 = %e\n",FEs,gbestfitness);
        end
        
        %% Birth of a new coyote from random parents (Eq. 7 and Alg. 1)
        parents         = randperm(n_coy_aux,2);
        prob1           = (1-Ps)/2;
        prob2           = prob1;
        pdr             = randperm(D);
        p1              = zeros(1,D);
        p2              = zeros(1,D);
        p1(pdr(1))      = 1; % Guarantee 1 charac. per individual
        p2(pdr(2))      = 1; % Guarantee 1 charac. per individual
        r               = rand(1,D-2);
        p1(pdr(3:end))  = r < prob1;
        p2(pdr(3:end))  = r > 1-prob2;
        
        % Eventual noise
        n  = ~(p1|p2);
        %生成子代小狼
        % Generate the pup considering intrinsic and extrinsic influence
        pup =   p1.*coyotes_aux(parents(1),:) + ...
            p2.*coyotes_aux(parents(2),:) + ...
            n.*(VarMin + rand(1,D).*(VarMax-VarMin));
        
        % Verify if the pup will survive
        pup_cost    = FOBJ(pup',fid);
        FEs=FEs+1;
        
        worst       = find(pup_cost<costs_aux==1);
        if ~isempty(worst)
            [~,older]               = sort(ages_aux(worst),'descend');
            which                   = worst(older);
            coyotes_aux(which(1),:) = pup;
            costs_aux(which(1),1)   = pup_cost;
            ages_aux(which(1),1)    = 0;
        end
        if gbestfitness> min(costs)
             gbestfitness= min(costs);
        end
        gbesthistory(FEs)=gbestfitness;
        fprintf("COA 第%d次评价，最佳适应度 = %e\n",FEs,gbestfitness);
        
        %% Update the pack information
        coyotes(packs(p,:),:) = coyotes_aux;
        costs(packs(p,:),:)   = costs_aux;
        ages(packs(p,:),1)    = ages_aux;
    end
    
    %% A coyote can leave a pack and enter in another pack (Eq. 4)
    if n_packs>1
        if rand < p_leave   %种群个体越多个体交换概率越大
            rp                  = randperm(n_packs,2);
            rc                  = randi(n_coy,1,2);
            aux                 = packs(rp(1),rc(1));
            packs(rp(1),rc(1))  = packs(rp(2),rc(2));
            packs(rp(2),rc(2))  = aux;
        end
    end
    
    %% Update coyotes ages
    ages = ages + 1;
    
    %% Output variables (best alpha coyote among all alphas)
    % [gbestfitness,ibest]   = min(costs);
    GlobalParams        = coyotes(ibest,:);
    %  gbesthistory(i) = gbestfitness;
    % fprintf('COA__%d---%e\n',i,gbestfitness);
    i=i+1;
    
    if FEs>=MaxFEs
        break;
    end
    
end
if FEs<MaxFEs
    gbesthistory(FEs+1:MaxFEs)=gbestfitness;
else
    if FEs>MaxFEs
        gbesthistory(MaxFEs+1:end)=[];
    end
end
end
